####Helper functions####

#Run BRGLM model for each topic separately
run_model_4_strong <- function(i){
  e1 <- myTryCatch(model <- brglm(off_rob~scale(get(i))+off_rob_l+factor(website_id)+
                                    factor(date_pos_rob),data=df_merged_mean,
                                  family=binomial(link="logit"),method="brglm.fit",pl=T))
  vcov <- vcovCL(model, cluster = ~ website_id,type="HC1")
  rob <- coeftest(model ,vcov = vcov) #Clustered SE
  se <- rob[,2] #Save SE
  p_value <- rob[,4] #Save P-value
  sign <- ifelse(p_value[2]<0.05/281,1,0) #Check whether p-value is smaller than Bonferreti adjusted p-value
  pos <- ifelse(rob[2,1]>0,1,0) #Check whether relationship is positive
  error <- ifelse(e1[[2]][1]=="Iteration limit reached",1,0) #Save whether there has been an error
  return(list(model,se,p_value,sign,pos,error))
}

#Simulations for marginal effects
marginal_model_sim <- function(variable_name){
  library(brglm) #Need to be loaded again inside due to parallel processing
  library(margins)
  library(sandwich)
  library(tidyverse)
  set.seed(123)
  formula <- as.formula(paste0("off_rob~",paste0(variable_name,"+off_rob_l+factor(website_id)+factor(date_pos_rob)")))
  m <-  brglm(formula,data=df_merged_mean,
              family=binomial(link="logit"),method="brglm.fit",pl=T)
  margins <- margins(m,vcov=vcovCL(m,cluster = ~ website_id,type="HC1"),vce="simulation", iterations=1000,
                     change="minmax",data=df_merged_mean)
  return(list(margins))
}

#Parallelizing for marginal effects models
run_marginal_models <- function(names){
  # Run marginal models in parallel
  # Calculate the number of cores
  no_cores <- detectCores(all.tests = FALSE, logical = TRUE)
  no_cores <- no_cores-1
  
  # Initiate cluster
  cl <- makeCluster(no_cores)
  
  clusterExport(cl, list("df_merged_mean","marginal_model_sim"))
  
  marginal_models <- parLapply(cl, 
                               c(names), 
                               marginal_model_sim)
  
  stopCluster(cl)
  return(marginal_models)
}


#Plot function for marginal effects
plot_simulations <- function(marginal_models_pos,
                             pos_names,names_output){
  marginal_models_eco  <- marginal_models_pos
  names <- pos_names
  
  
  for (i in seq(names)){
    
    name <- names[i]
    names_output_temp <- filter(names_output,names_fe %in% name)
    #Get position of variable name
    pos1 <- which(summary(marginal_models_eco[[i]][[1]])[,1] %in% name)
    marginal_frame_pos <- data.frame(Variable = names_output_temp$topic_names_text_fe,
                                     Category = names_output_temp$category,
                                     Coefficient = summary(marginal_models_eco[[i]][[1]])[pos1,2],
                                     SE =  summary(marginal_models_eco[[i]][[1]])[pos1,3],
                                     Model = "FE") 
    if (i==1){
      marginal_frame <- marginal_frame_pos
    }
    else {
      marginal_frame <- rbind(marginal_frame,marginal_frame_pos)
    }
    rm(marginal_frame_pos)
  }
  
  
  #Confidence intervals
  interval1 <- -qnorm((1-0.95)/2)  # 95% multiplier, non-Bonferetti-corrected
  interval2 <- -qnorm((1-(1-0.05/281))/2)  # 95% multiplier, Bonferetti-corrected
  
  plot_pos <- ggplot(marginal_frame, aes(colour = Model)) +  scale_colour_grey(start = 0, end = 0.8) 
  plot_pos <- plot_pos + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
  plot_pos <- plot_pos +  geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                             ymax = Coefficient + SE*interval1),
                                         lwd = 1.5, position = position_dodge(width = 0.9))
  plot_pos <- plot_pos + geom_linerange(aes(x = Variable,  ymin = Coefficient - SE*interval2,
                                            ymax = Coefficient + SE*interval2),
                                        lwd = 1/2, position = position_dodge(width = 0.9))
  plot_pos <- plot_pos + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                             ymax = Coefficient + SE*interval2, shape = Model),
                                         lwd = 1/2, position = position_dodge(width = 0.9)) + scale_shape(solid = TRUE)
  plot_pos <- plot_pos + coord_flip() + theme_bw() +theme(axis.text=element_text(size=14,face="bold"),
                                                          axis.title=element_text(size=14,face="bold")) 
  plot_pos <- plot_pos + ggtitle("") + xlab("") + ylab("") + scale_size(guide = 'none') + 
    facet_wrap(~ Category, nrow=3, scales = "free_y") +
    theme(legend.text=element_text(size=12,face="bold"),
          strip.text = element_text(size=12),
          legend.position="none",
          panel.spacing.x = unit(10, "mm"))
  return(plot_pos)
}


#Define function to catch error messages 
myTryCatch <- function(expr) {
  warn <- err <- NULL
  value <- withCallingHandlers(
    tryCatch(expr, error=function(e) {
      err <<- e
      NULL
    }), warning=function(w) {
      warn <<- w
      invokeRestart("muffleWarning")
    })
  list(value=value, warning=warn, error=err)
}